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ABSTRACT 

In this paper, a quantitative characterization for the evolutionary sequence of stel- 
lar self-gravitating system is investigated, focusing on the pre-collapse stage of the 
long-term dynamical evolution. In particular, we consider the quasi-equilibrium be- 
haviors of the A^-body systems in the setup of the so-called Antonov problem, i.e., 
self-gravitating A-body system confined in an adiabatic wall and try to seek a pos- 
sible connection with thermostatistics of self-gravitating systems. For this purpose, a 
series of long-term A^-body simulations with various initial conditions are performed. 
We found that a quasi-equilibrium sequence away from the thermal equilibrium can be 
characterized by the one-parameter family of the stellar models. Especially, the stellar 
polytropic distribution satisfying the effective equation of state P oc p 1+1 / n provides 
an excellent approximation to the evolutionary sequence of the A^-body system. Based 
on the numerical results, we discuss a link between the quasi-equilibrium state and 
the generalized thermostatistics by means of the non-extensive entropy. 
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^ ; 1 INTRODUCTION 



1.1 Self-gravitating A-body systems and thermostatistics 

The long-term dynamical evolution of stellar self-gravitating system driven by the two-body relaxation is an old problem with 
rich history in astronomy and astrophysics and even in statistical physics. The problem, in nature, involves the long-range 
attractive nature of gravity and because of its complexity and peculiarity as well as the physical reality, astronomers and 
statistical physicists have attracted much attention on this subject. 

Historically, an important consequence from the thermodynamical arguments had arisen in the 1960s. lAntonovl dl962t) 
first pointed out that no stable equilibrium state exists for a high-dense clusters. To prove this, he considered a very idealized 
situation called Antonov problem, i.e., a stellar self-gravitating system confined in a spherical cavity with radius r e (see Fig0. 
Then, under keeping the energy E and the mass M fixed, the standard statistical mechanical approach based on the maximum 
entropy principle leads to the c onclusion that no stable e quilibrium state exists for the larger radius r e > A cr i t (-E/GM 2 ), 
where A cr it is the critical value iPadmanabhaniri989l Il990t for pedagogical reviews) . Note that the numerical value of A cr i t 
depends on the choice of the entropy and Antonov obtained A cr i t = 0.335 in the case adopting the Boltzmann-Gibbs entropy: 



Sbg = — 



d s xd 3 v f(x,v) In f(x,v), 



(1) 



where f(x,v) is the one-particle distribution function in phase-space. Later, iLvnden-Bell fc Woodl ( 1968|) re-examined this 
issue and showed that the unstable thermal state found by Antonov can be explained by the thermodynamic instability arising 
from the negative specific heat. They especially called the instability gravothermal catastrophe. 

Since the 1960s, the gravothermal instability discovered by Antonov has become a standard notion of the stellar dynam- 
ics and the role of the instability has been extensively discussed. Thanks to the unprecedented development of the computer 
facility as well as the sophisticated numerical techniques based on the A-body simulation and the Fokker-Planck calcula- 
tion, our view of the late-time phase of the stellar gravitating system has dramatically improved (e.g.. lMevlan fc Heggidll997l : 
iHeggie fc Hutl2003|) . Among various theoretical de velopments, one important lan dmark wo uld be a discovery of the qravo ther- 
mal oscillation, which was originally suggeste d bvlSugimoto fc Bettwieserl (119831) (see also iBettwieser fc Sugimotclll984|) and 
was later confirmed by numerical simulation jMakincE)96l) ""with a great advantage of a special purpose hardware, GRAPE 
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(e.g., ISugimoto et~al"lll99Cl ; iHut fc Makinolll999li . the core-collapse triggered by gravothermal instability was shown to be 
terminated by the formation of binary as a result of three-body interaction and the oscillatory behaviors of the core motion 
has been clearly revealed. The gravothermal oscillation is thus thought to be one of the most fundamental stellar dynamical 
processes and provides a basis to understand the dynamical history of globular clusters as real astronomical system. 

So far the theory of long-term evolution in stellar self-gravitating syst ems has been dev elop ed without any recours e to th e 
statistical mechanics or thermodynamics except for the seminal works bv lAntonovl il962t) and iLvnden- Bell fc Woo J £L968). 
This is essentially because the issues under consideration are, in nature, non- equilibrium prob lem with long-range interaction 
and usual sense of the thermal equilibrium becomes inadequate. Indeed, as lAntonovl dl962l) emphasized, no strict meaning 
of the thermal equilibrium exists in the case of an isolated stellar system without boundary. Nevertheless, thermostatistical 
point-of-view is helpful in predicting the fate of the unstable system. This would be even true in the non-equilibrium situation 
as long as the collisional evolution during the relaxation timescales is concerned. In this sense, there might be some possibilities 
to recast the non-equilibrium self-gravitating systems in terms of an extended view of the thermostatistics. The present paper 
attempts to address these issues by considering an old problem of stellar dynamics from a somewhat new point-of-view. 
Especially, we wish to focus on the thermostatistical characterization of non-equilibrium self-gravitating system by means of 
the new thermostatistical formalism and discuss the validity of it. 

Indeed, the self-gravitating system may be the best testing grou nd for the rece nt postulated introduction of a non-extensive 
generalization of Boltzmann-Gibbs statistics, originally proposed by Tsallis ( 1988). In contrast to the Boltzmann-Gibbs entropy 
Q which is properly an extensive variable, the entropy used in the 'Tsallis' formalism is pseudo-additive and/or non-extensive 
one 1 : 

S q = -^-j J d 3 xd 3 v [{p(x, v)} q - p(x, v)] , (2) 

where the probability p(x,v) is not one-particle distribution function but a sort of primitive phase-space distribution that 
satisfies the normalization condition: 

Note that the Boltzmann-Gibbs entropy Q is recovered in the limit q — > 1 and in this limit, the primitive distribution is 
identified with the one-particle distribution. 

Because of its non-extensive nature, the non-extensive 'Tsallis' formalism might have a potential to deal with the long- 
range systems including the self-gravitating systems. Further, it is expected to deal with a variety of interesting non-equilibrium 
problems such as quasi-steady state or quasi-cquilib rium state far from the the rmal equilibrium, to which the standard 
Boltzmann-Gibbs statistics cannot be applied jTsallislll999l: l Abe fc Oka moto 200lJ, for comprehensive reviews). Nevertheless, 
most of the works on this subject are concerned with construction of a consistent formal framework and little works have been 
known concerning the physical realization of the non-extensive statistics. Hence, it seems interesting to address the issues on 
the reality of non-extensive statistics. 



1.2 Non-extensive generalization of the thermostatistical description of self-gravitating system 

Among various issues on the physical rea lity of the non-extensive thermostatistics, we have recently re-examined the clas- 
sic problem considered by lAntonovl l|l96 ^) bv means of the non- extensive thermostatistics with Tsallis' generalized entropy 
jTaruva fc Sakaeamil l200l l2003allbld: ISakaeami fc Taruvall2004 for a review). According to the framework using normal- 
ized q-expectation values bv lT^alii^llMel'Kki^^Pl^^ the one-particle distribution f( x,v) is expressed in terms 
of the primitive distribut ion p{x,v) by (For the cases using standard linear-mean values, see IPlastino fc Plastind Il993l : 
iTaruva fc Sakagamj|2002T) : 

/(«,«) = M r . (4) 

d 3 xd 3 v {p(x,v)} q 

Then, the standard analysis of statistical mechan ics based on the Tsallis ent ropy @ reveals that the extremum state of the 
entropy is reduced to the power-law distribution (ITaruva fc SakagamiE)03bT) : 



f(x,v) =A 



1 2 

$0 - -v - $(x) 



9/(1-9) 



(•») 



where A and $0 are the numerical constants and $(a;) is the gravitational potential. In terms of the specific energy of a 
particle defined by e = v 2 /2 + $(x), this gives / oc ($0 — e) q ^ 1 ~ q \ Equation <|5)l is the famo us one-parameter family of the 
stellar models referred to as the stellar polytropic distributio n (e.g.. iBinnev fc Tremainelll98^ and the g-parameter is related 
to the polytrope index given by iTaruva fc Sakagamil l2003b'l : 

1 Strictly speaking, equation is not exactly the same entropy as originally introduced by Tsallis in a sense that the entropy is 
defined in /i-space, not in T-space. Nevertheless, we keep to call it Tsallis entropy bec ause the thermostatistical f ramework with that 
is constructed self-consistently is quite similar to the standard Tsallis formalism (see ITaruva fc Sakagami 200 3bh . 
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Setup of the Antonov problem. 
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Similar to the case of the Boltzmann-Gibbs entropy (0, the statistical mechanical analysis based on the non-extensive entropy 
also showed that there exists the thermodynamic instability, which can be consistently explained from the presence of 
the negative specific heat in terms of the non-extensive thermodynamics. Thus, the stellar polytropes as extremum states 
of Tsallis entropy are consistently characterized by the non-extensive thermostatistics and may play a role of the thermal 
equilibrium and/or quasi-equilibrium like the extremum state of the Boltzmann-Gibbs entropy. Nevertheless, concerning their 
reality, it is still unclear whether the consistent thermostatistical results really imply the existence of quasi-equilibrium state 
or not. In order to get a further insight into this issue, apart from the thermostatistical analysis, dynamical and/or kinematical 
aspects of the self-gravitating system should be investigated in details using the A -body simulation, which we will discuss 
here. Some of the work we report in this paper has already appeared in Letter form (iTaruva fc SakagamjE)03cl) . This paper 
further addresses new additional JV-body experiments and an improved analysis, together with some careful discussions. 



1.3 Outline of this paper 

In section [5] within a setup of the classic issue considered bv lAntonovl Jl962t) and lLvnden- Bell fc WoodlJl96^ . we first address 
some important properties of the quasi-equilibrium distribution characterized by non-extensive entropy by comparing with 
those of the standard Boltzmann-Gibbs statistics. We then move to discuss the A-body simulations. In section|21 the A-body 
treatment of the Antonov problem is considered and the initial conditions are summarized. Section [I] is a main part of this 
paper, which describes the results of A^-body experiments in details. After checking the A'-body code in section |T?T1 we discuss 
the quasi-equilibrium behaviors of the long-term evolution starting with the stellar polytropic distribution in section FOl and 
try to fit their evolutionary sequences by the stellar polytropes, which are shown to be a remarkably good fit. The results are 
also compared with an alternative one-parameter family of stellar model, i.e., the King model. In section FOl we attempt to 
clarify the condition for quasi-equilibrium state characterized by the stellar polytropes in other class of initial conditions, i.e., 
the stellar models with cusped density profile. Finally, section^] is devoted to the discussion and the conclusion. 



2 ANTONOV PROBLEM AND ITS NON-EXTENSIVE GENERALIZATION 

In this section, we briefly mention the Antonov problem and the basic properties of the equilibrium states characterized by 
the extensive and the non-extensive entropies. 

First recall the setup of the so-called Antonov problem. In this problem, we consider the many-body particle system 
confined in an adiabatic wall (see Fig0. The radius of the wall is given by r e and we assume that the total mass M and the 
total energy E are kept fixed. All the particles in this system interact via Newton gravity and bounce elastically from the 
wall. For further simplification, each particle is assumed to have the same mass m. Under these conditions, we consider a fully 
relaxed system and clarify the nature of stability or instability for such a system. Note that there are two characteristic time- 
scales in self-gravitating system. One is the dynamical time tdyn and the other is the relaxation time treiax- Usually, for a system 
conta ining a large number of particles N 1, these time-scales scale as t rc iax — (0.1AT/ In N) td yn (e.g., Bin nev fc Tremaind 
119871) . Accordingly, the equilibrium states considered here should be a long-lived state much longer than t re i ax and thus the 
thermostatistical consideration becomes helpful. 

In a language of thermostatistics, a fully relaxed equilibrium state corresponds to the state that maximizes the entropy. 
If we adopt the standard Boltzmann-Gibbs entropy the equilibrium distribution is obtained by extremizing the entropy 
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Sbg under keeping the mass M and the energy E fixed: 

M = / d 3 xd 3 v f(x,v), (7) 
E = 

Then t he first variation of entropy with respect to /, i.e., ^[6*bg — a M — f3 E] = with a and (3 being Lagrange multipliers 
~lll9fi " ' 

3/2 



1 2 1 _ j \ ' 

-v 2 + -$>{x) 



fi.v.v) : 4>(,ri = -G / dVdV {^T^TT- (8) 



then the hrst variation ot entropy with respect to /, i.e., o|obg — a iw 
yields ijAntonovlll'ffil iLvnden-Bell fc Woo Jl968UPadmanabhadll98lih : 

/(*,») = (^) 3/2 ^) c_/3 " 2/2 ( 9 ) 
with p being the mass density of the system defined by 

p(x) = J d 3 v f(x,v). (10) 

Computing the pressure, one can deduce that equation of state is indeed isothermal in the case of the isotropic velocity 
distribution: 

P(x) = J d A v ±v 2 f(x,v), (11) 

which leads to P(x) = P~ 1 p(x). 

On the other hand, if one uses the non-extensive entropy the variational problem with respect to the primitive 
distribution p(x,v) not the one-particle distribution f(x,v) finally leads to the stellar polytropic distribution ©. With a 
help of the relation @, this distribution can be rewritten with 

f(xv) = i eM u m r 3/2 (12) 

n,} 4^5(3/2,71-1/2) {(n+l)P(x)/p(x)}s/a \ ( n + 1) P(x)/p(x) J ' V ' 

where the function B(a, b) represents the beta function. Note that we used the following relations to derive the above expression: 

p(r) = 4V2tyB Y~~^j A ~ *( r )] 1/<1 ~ 9)+1/2 ' (13) 

P(r) = ^flfiAVl^-^p-^ 2 (14) 



v 2 1 

with r being the radius r = \x\. These two equations lead to the polytropic equation of state: 

P(r) = K n p 1+1/n (r). (15) 

The explicit form of the dimensional constant K n is given in iTaruva fc Sakagami J2003bll (Eq.(21) of their paper). 

Comparing the distribution function with ©, the stellar polytropic distribution can be reduced to the isothermal 
distribution when taking the limit n — > +oo, and only in this limit, the isothermal equation of state is recovered. Basically, 
the equilibrium properties of the stellar polytropes are similar to those of the isothermal distribution. But, focusing on their 
thermodynamical stabilities, there exists some distinctive features. Before showing this, we first note that the one-particle 
distribution function 1121 does not yet completely specify the equilibrium configuration, because functional forms of the 
density and the pressure are not determined. To determine the equilibrium configuration, we need to specify both p(r) and 
P(r), or equivalently the gravitational potential $(r). This can be achieved by solving the Poisson equation explicitly. Under 
the spherically symmetric configuration, the Poisson equation becomes 

i d f r idm\ =4irGp{r) (16) 



r 2 dr 1 dr 
or the condition of hydrostatic equilibrium: 

dP(r) Gm(r) . dm(r) , . 2 /,„n 

-fiT = ^PM> - ^P(r)r, (17) 

where the quantity m(r) denotes the mass inside the radius r. Denoting the central density and pressure by p c and P c , we 
then introduce the dimensionless quantities: 

1/2 



P = Po[0(O] n , r = <j ^ - ; } C, (18) 



(n + l)P e 
AnGpl 

which yields the following ordinary differential equation 



+ ^6' + 6 n =0, (19) 
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where prime denotes the derivative with respect to £. To obtain the physically relevant solution of 11911 . we put the following 
boundary condition: 

0(0) = 1, 0'(O) = 0. (20) 

A family of solutions sat isfying 1200 is the Emden solution, which is well-known in the subject of stellar structure (e.g., see 
chapter IV o f 1Chandrasekharll93aT To characterize the equilibrium properties of Emd en solutions, it is convenient to introduce 
the f ollowing set of variables, referred to as homology invariant variables (e.g.. Ichandrasekhajll939l: iKippenhahn fc Wigertl 
Il990ll : 

rflnm(r) _ 4?rr 3 p(r) _ £6 n 

U = dlnr ~ m(r) ~ __ r ' ( } 

dlnr P{r) r y ' 9 ' y ' 

which reduce the degree of equation I19H from two to one. One can evaluate the total energy of the confined stellar system in 
terms of the pressure P e ,the density p e at the boundary r e and the total mass M: 

„ 3 f r " , 2 „, ■. f ,e , Gm(r) dm 

2 J **^P<r)-J o * 



i 



3 f GM , , N MP e 
2 



l)HLLi K + ( n -2)4ivr 3 e P t 



by which the dimensionless quan tity A can be expressed as a function of the homology invariant variables at the wall 
jTaruva fc SakagamilEo02l.l2003allbl) : 



A: Er * 



2 I V e J V 



In figure [5] the dimensionless quantity A is plotted as a function of the density contrast D = p c / p e , i.e., the ratio of 
the central density to that at the boundary. Each curves in (A, D)-plane corresponds to the equilibrium sequence of stellar 
polytrope with a different value of polytrope index. For comparison, we also plot the equilibrium sequence of isothermal 
distribution in thick solid line. Figure [5] shows that the equilibrium sequences can be classified as two types. One is the the 
equilibrium sequences with indices n ^ 5. In this case, the A-curves monotonically increase with density contrast. On the 
other hand, the other type of equilibrium sequences with n > 5 have peaks in each trajectory. At a certain value of the density 
contrast D C rit, the trajectory reaches a maximum A cr i t and beyond this, the monot onicity is lost. This is even true in the limit 
n — > +00 and the m aximum value of A and the critical value of D can be read off iAntonovll962l : iLvnden- Bell fc WooJll96gt 
|Padmanabhanlfl98St) : 

A crit = 0.335, Dcrit = 709. (24) 

In the case of Boltzmann-Gibbs entropy, it had been already known that, along the curve X(D) derived from the condition 
<5<Sbg = 0, all states with D > -D cr it are unstable. This can be proved by th e turning-point analysis for a linear series of 
equilibria (iLvnden-Bell fc Woodlll96l |Padmanabhanlll989t lKatjll978t 1 19791) . Al so, the same results can b e obtained by 
explicitly evaluating the eigenmodes for the second variation of the entropy 8 2 Sb.g JPadmanabhar]ll990l. Il989^ . The absence 
of stable thermal equilibria in this regime clearly indicates the instability, referred to as the the gravothermal catastrophe. 

In similar manner to the isothermal distribution, it has been recently shown that the stellar polytrope confined in an 
adiabatic wall exhibits the gravothermal instability in the case of the poly tropic index n > 5. The evaluation of eigenvalues 
for the second variati on of the Tsallis entropy S 2 S a a lso derives the same results as the above turning-point analysis in 
terms of the A-curvc jTaruva fc Sakagamil 12003 . l2003allbl) . Note that the gravothermal instability in the isothermal case is 
heuristically explained by the presence of negative specific heat. That is, the co-existence of the regions with negative and 
positive specific heat leads to the violation of thermal balance and finally causes the catastrophe temperature growth. This 
can be achieved when the radius of the adiabatic wall or equivalently the density contrast becomes sufficiently larger so that 
the effect of self-gravity is important at the inner dense region and becomes negligible at the outer sparse region. Although 
this intuitive explanation heavily relies on the standard notion of extensive t hermostatistics, it is empha sized that the same 
interpretation can be possible in the case of the stellar polytropic distributions (iTaruva fc Sakagami2003ah . Therefore, despite 
some distinctive features, the stellar polytropes are consistently characterized as equilibrium distribution by means of the non- 
extensive formalism and their thermodynamic properties are similar to those of isothermal distribution. 

Finally, the above arguments indicate that the stellar polytropic distribution might be regarded as a sort of thermal 
equilibrium, since it is the extrema of the Tsallis entropy. However, polytropic equation of state for the stellar polytropes 
clearly shows that the one-dimensional velocity dispersion defined by 

p{r) 

manifestly depends on the radius r (see Eg. 1151 1. Only in the isothermal case n — > oo, (j Vj id is kept spatially constant. Thus, 
it is expected that a gradient of the velocity dispersion is relaxed on timescales of two-body encounter. This means that the 
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D(-P> e ) 

Figure 2. Energy-density contrast relationship of stellar polytropes and evolutionary tracks of the quasi-equilibrium states observed in 
the TV-body experiments. Each curve represents equilibrium sequence of stellar polytrope with a different polytrope index n. The thin 
dotted line which traces the turning points for each equilibrium sequence is the marginal stability boundary inferred from the Tsallis 
entropy S q . The thick arrows show the evolutionary tracks of TV-body simulation starting with stellar polytropic and/or isothermal 
distributions (see table HI and l2l, 

stellar polytrope is no longer the equilibrium but quasi-equilibrium state. At present, it might be also possible to give another 
interpretation that the stellar polytropic states with n < 5 an d with n > 5 and D < D cr it are dyna mical equilibrium states 
that are even stable with respect to a non-linear perturbation llChavanisll2003t Ichavanis fc Sirell2004f) . In subsequent section, 
we report the results of the TV-body simulations, which are carried out to investigate how the stellar polytrope actually evolves. 



3 NUMERICAL SIMULATIONS 

We are in a position to discuss the long-term dynamical evolution from the non-equilibrium TV-body system and to explore 
the reality of stellar polytropes as quasi-equilibrium state. The TV-body experiment considered here is the same situation as 
investigated in classic papers. That is, we treat the self-gravitating TV-body system confined in an adiabatic wall of sphere. 
Hereafter, we set the units to G — M = r e = 1 without loss of generality. According to the Antonov problem, all the particles 
are assumed to have the same mass m = 1/TV, where TV is the total number of the particles. The initial conditions for particle 
data were whole created by a random realization of the stellar models. To investigate the evolutionary states of the TV-body 
system from the thermostatistical point-of-view, we deal with the three kinds of the stellar models: (i) isothermal distribution 
(ii) stellar polytropes and (iii) the stellar models with cuspy density profile. Table [5] and [3] summarize the simulation 
runs, which are frequently referre d in subsequen t sectio n. The generation o f random initi al data was basically followed by 
the rejection method described in lAarseth et all <ll974l) (see also Chap. 8 of lAarseth 2003). Using the analytic form of the 
mass profile m(r) and the one-particle distribution function /(e) as function of specific energy e = v 2 /2 + 4>(r), this method 
generates the random distribution that possesses the same phase space structure as in the initial stellar model. 

We are specifically concerned with collisional aspects of TV-body dynamics, the timescale of which is much longer than the 
two-body relaxation time. For this purpose, we utilized a special-purpose har dware, GRAPE-6, w hich is especially designed to 
accelerate the gravitational force calculations for collisional TV-body system iMakino et alJl^OOot) , With this implementation, 
we use the fourth-order Hermite integration scheme with individual time-step algorithm, which is suitably efficient to combine 
with GRAPE-6 facility. 

In our setup, the adiabatic wall is implemented by the same procedure as used bv lEndoh et alJ |l997). The adiabatic 
wall reverses the radial components of the velocity for particles just located at the wall. Since we use an individual time-step 
algorithm, the positions of the particles are monitored at each fixed time interval AT, which is chosen to be a synchronized 
time interval in our time-step algorithm. This guarantees that the penetration of particles into the wall can be ignored. At 
these times, the radial components of velocities for particles outside the wall are reversed if the radial velocity vector is directed 
outward. Namely, 
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(v ■ x) X 

v — > v — —. — - — - (26) 

\x\ \x\ 

After particles were reversed, we recalculated the force and adjusted the time-step of particles. 

Note that we did not use the regularization scheme to treat the close two-body or multiple encounter f e.g. . lAarsetbl2003f) . 
Because of the adiabatic wall, the standard scheme of regularization method is not directly applicable and the implementation 
of the regularization keeping the high accuracy requires a considerable amount of programming. Since we are interested in a 
quasi-equilibrium evolution before the core-collapse stage, absence of the regularization scheme itself is not crucial. Rather, 
a serious problem might arise from an introduction of the potential softening in order to reduce the numerical error. In 
general, the potential softening diminishes the close-encounter of each particle and it makes the energy exchange of the 
particles inefficient. This would lead to the overestimation of the relaxation timescales even before the core-collapse stage. In 
appendix A, the significance of the potential softening to the time-scale of core-collapse and/or quasi-equilibrium sequences 
is examined. Based on these experiments, we adopt the Plummer softened potential (<f> = l/\/r 2 + e 2 ) and basically set the 
softening parameter to e= 1/N when estimating the time-scale (Sec l4.2ll . Otherwise, we set e = 4/iV (Sec l4.1l and 14. 31 . 2 



4 RESULTS 

In our present situation with units G = M = r e = 1, the dynamical time given by roughly corresponds to t^yn = 
(37r/16Gpo) 1//2 — 1 and thus the global relaxation time becomes t re iax = 0.1(N/ InN) tdyn ~ N/\nN. While this relation 
gives a crude estimate of the relaxatio n timescales, a more useful convention might be the half-mass relaxation time (e.g., 
ISpitzedir987l: iBinnev fc Tremaind[t987h : 

N / r 3 V /2 

«* = °- 138 inA [mi ) • (27) 

where the Coulomb logarithm InA is usually taken as InA = ln(0.4/V) or In(O.LZV) ljSr>itzerlll987l : lOiersz fc HeggieHl994|) . In 
what follows, adopting the latter convention (A = 0.1JV), all the simulation results are presented by rescaling the timescale 
with the half-mass relaxation time evaluated at an initial time, t r h,i. In tables Hl3l half-mass radii for each initial distribution 
are evaluated and their numerical values are summarized. 



4.1 Isothermal distribution 

Let us first check our numerica l calculations by exa mining the cases with isothermal initial conditions and comparing those 
results with previous work by lEndoh et alJ <ll997T) , who have investigated the gravothermal expansion of the isothermal 
distribution under the same setup as examined in our simulation. Following their paper, we examine the three kinds of initial 
conditions, summarized in table 

Figure|H]andg]respectively show the snapshots of evolved density profiles and one-dimensional velocity dispersion profiles 
starting from a stable initial configuration (left) and unstable initial distributions(midrf/e, right). Here, the term, stable or 
unstable is used according to the initial density contrast D smaller or larger than the critical value D C rit = 709, which is 
obtained from the thermostatistical prediction. The evolutionary path for each run is depicted in figure |5] In figure |HJ while 
the evolved density profiles starting from the the stable initial condition is, by construction, stable and almost remain the 
same, the fate of the unstable cases sensitively depend on the randomness of the initial data, leading to the very different 
endpoints. The results depicted in the middle and the right panels of figure^] both of which were started with the same initial 
parameters, are indeed such examples. Discriminating the final fate of the system from the initial random distributions seems 
generally difficult, however, the evolved density profiles seen in figure [3] are intimately connected with the response of the 
velocity structure shown in figure 0] Suppose that the velocity dispersion at the central part is initially higher than that at 
the outer part. In this case, the energy exchange by the two-body encounter yields the kinetic energy transport from the core 
to the halo. In other words, the outward heat flows occur and due to the negative specific heat, the inner part gets hotter, 
leading to the catastrophic growth of velocity dispersion. As a consequence, the core-collapse takes place. On the other hand, 
when the velocity dispersion at the inner part is relatively lower than that at the outer part, the inward heat flow conversely 
occur, which makes the core expand. In contrast to the former case, the heat flow can stop after balancing the thermal inertia 
at inner and outer parts. As a result, the system finally reaches the stable isothermal configuration. Though not clearly seen 
its signature, the response of the velocity dispersion profiles seen in figure [I] are c onsistent with the e volution of the density 
profiles. Compared the results in the unstable cases with the results obtained bv lEndoh et alJ dl997f) . the resultant density 
profiles are similar and the evolved timescales for the unstable cases almost agree with each other. 

As for the stable case, a calculation was proceeded up to t = 112£ r h,i, corresponding to t = 2000 in unit of A-body time, 
which was longer calculation time than that of the other runs. A closer look at figure [3] shows that the central density slightly 
increases and the system seems to become out of equilibrium. The total energy of the system was conserved to 0.05% at the 
end of the calculation, which is relatively larger value than the errors in the other runs. Thus the destabilization of the core 

2 Hereafter, the quantity e is interchangeably used to imply the specific energy of a particle and the softening length. Readers should 
not confuse the meanings of e. 
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Figure 3. Snapshots of density profile in the stable case(Le/i) and the unstable ca,ses(Middle and Right). In each panel, the thick lines 
represent the snapshots of simulation results, while the thin lines indicate the ones obtained from the Emden solutions of isothermal 
distribution. 
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Figure 4. Snapshots of one-dimensional velocity dispersion profile in the stable case(Left) and the unstable cases(Middle and Right). 



Table 1. Model parameters of initial condition in the case of the isothermal distribution 



run # 


parameters 


half-mass radius (rh/r e ) 


# of particles 


status 


stable 


D = 110 


0.4834 


2k 


stable 


unstablel 


D = 10 4 


0.4993 


2k 


expansion 


unstable2 


D = 10 4 


0.4993 


2k 


collapse 



may be attributed to the numerical error in the energy conservation. Though the effect of destabilization yields a serious 
problem in estimating the time scale of core-collapse, we are specifically concerned with the evolutionary sequence before 
entering the core-collapse phase. Hence, at a level of this present accuracy, one can ignore the destabilization effect as long as 
the total energy of the system is well conserved. 



4.2 Stellar polytropes 

Having checked the validity of the numerical calculations, we next investigate the non-equilibrium evolution away from 
the isothermal equilibrium state. In this subsection, we deal with a class of initial conditions characterized by the stellar 
polytropes. Table [5] summarizes the model parameters of initial distributions. The evolutionary track for each run is depicted 
in the energy-density contrast relation in figure 

Overall behaviors of the simulation results are as follows. Before entering the core-collapse stage, most of the system 
experiences the quasi-equilibrium stage, in which the distribution function slowly changes in time. The evolutionary sequence in 
the quasi-equilibrium stage can be well-approximated by the one-parameter family of stellar polytropes with the time- varying 
polytrope index. Apart from some fluctuations, the fitted value of the polytrope index, on average, increases monotonically 
with time, implying that the system tends to approach the exponential distribution. 

To see the quasi-equilibrium behavior quantitatively, we first show the representative result obtained from the run n3A, 
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Table 2. Model parameters and the evolutionary states in the cases starting from the stellar polytropes 



9 



run # model parameters half-mass radius (rh/r e ) # of particles realizations fitting to stellar polytrope 



n3A 4 


n = 


3: 


D = 


10 4 


0.3282 


2k, 4k, 8k 


2, 1 ,1 


successful 


n5A 


n = 


5, 


D = 


10 2 


0.4611 


2k 


2 


successful 


n5B 


n = 


5, 


D = 


10 3 


0.3114 


2k 


3 


successful 


n5C 


n = 


5, 


D = 


10 4 


0.2025 


2k 


2 


successful 


n5D 


n = 


5. 


D = 


10 6 


0.08196 


2k 


2 


failed 


n6A 2 


n = 


6, 


D = 


110 


0.4531 


2k 


1 


successful 


n6B 3 


n = 


G, 


D = 


10 4 


0.1794 


2k 


1 


successful 



1 initial distribution corresponding to run A in|Tanrv^^^Sakaeami|( 2003c 

2 initial distribution corresponding to run Bl in Taruya & Sakagami ( 2003i 
initial distribution corresponding to run B2 in lTaruva & Sakagami ( 2003. 



i.e., the stellar polytrope with index n — 3 and with the initial density contrast D = 10 4 . Then we discuss the other runs, 
n5A ~ n5D and n6A, B in section 14. 2. 21 

4-2.1 runnSA 

Figure shows the time evolution of the Lagrange radii taken from the run with N = 2K, plotted as function of time in unit 
of the half-mass relaxation time. With the softening parameter e = 1/N, the core-collapse takes place at t ~ 44t r h,i and the 
core-halo structure was developed at the end of the calculation. Looking at an earlier phase, the Lagrangian radii evolve very 
slowly and one can clearly distinguish the timescales between the early relaxation phase and the late-time core-collapse phase. 
Since the early stage of the time evolution seems nearly equilibrium, we will especially call it quasi-equilibrium evolution. 
In a quasi-equilibrium regime, while decreasing the inner Lagrangian radii, the outer Lagrangian radii slightly expands to 
compensate the slow contraction of the core. 

To see the quasi-equilibrium structure in more detail, in figure [U we plot the snapshots of the density profile as function 
of radius (left), the distribution function as function of specific energy of the particles e = « 2 /2 + &(r) (middle) and the 
velocity dispersion profile as function of radius (right) during the quasi-equilibrium regime. In each panel, the symbols denote 
the simulation results. Note that for clarify, the results are offset vertically, successively by two-digits below, except for the 
final output at t = 30i r h,i- 

In figure |S| the solid lines show the initial stellar polytrope with n = 3 evaluated from the Emden solutions. Comparing 
those curves with simulation results, one deduces that the system at an earlier time is slightly out of equilibrium and it 
gradually deviates from the initial state. While the velocity dispersion profile monotonically increases at an outer part, the 
density profile first increases at both the inner and the outer parts (t ,$ 5f r h,i), leading to a slight decrease of the density 
contrast D = p c /p e (see the arrow labeled by n3A in figure |5J. Later, the increase of the core density surpasses that of the 
edge density and thereby the density contrast turns to increase (t £ 10t r h.i). 

In figure |S| in order to characterize the evolutionary sequence of the quasi-equilibrium state, the simulation results are 
compared with a sequence of the stellar polytropes. The fitting results are then plotted as long-dashed, short-dashed, dot- 
dashed and dotted lines from the data at t — 5t r i,,i to that at t — 30t r h,i- In fitting the simulation data to the stellar polytropes, 
we first quantify the radial density profile p(r) from each snapshot data. Selecting the 100 points from it at regular intervals 
in logarithmic scale of radius r, the results are then compared with the Emden solutions. 3 Note that in our present situation, 
the total energy and the mass are conserved. Thus, the only fitting parameter is the polytrope index. We determine the index 
n so as to minimize the function x 2 : A 

*(„\ - / Psim(rQ - PEmdcn(n; w) 1 ^ 

Clearly from figure |S| the stellar polytropes quantitatively characterize the evolutionary sequence of the simulation 
results. Note that we obtained ~ 4-8 for each time-step, indicating that the fitting result is satisfactory. A closer look 
at the low-energy part of the distribution function /(e) reveals that the simulation results partly resemble the exponential 
form rather than the power-law function. Nevertheless, the most remarkable fact is that the stellar polytropes as simple 
power-law distribution globally approximate the simulation results in a quite good accuracy. Moreover, the fitting results in 

3 Reason why we adopted the density profile instead of the distribution function is that the density profile can be very sensitive to 
the dynamical evolution process. Further, in spherically s ymmetric and isotropic cas es, the density profile is uniquely determined from 
the distribution function through the Eddington formula jBinnev fc Tre mainc 1987). Thus, in the present situation, it seems better to 
perform the fitting by using the radial density profile. 

4 Strictly speaking, this is not the x 2 function usually used in the likelihood analysis. 
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figure|S|indicate that the fitted values of the polytrope index monotonically increase in time, in contrast to the non-monotonic 
behavior of the density contrast. In order to quantify the evolution of polytrope index, we collect the new snapshot data at 
each time interval At = 10 in TV-body units for the run with TV = 2K. Repeating the same fitting procedure as in figure E] 
the polytrope indices are estimated at each time step and the resultant values are plotted in figure [7| together with the \ 2 
value of the fitting results. Here we also plot the results obtained from the run with TV = and TV = 8K, in which the time 
intervals are respectively chosen as At = 20 and 50. 

Clearly, the fitted values of the polytrope index monotonically increase apart from the fluctuations during the short time 
interval. The growth rates of the polytrope index normalized by half-mass relaxation time almost coincide with each other, 
consistent with the fact that the quasi-equilibrium sequence is evolved via two-body relaxation. In upper panel of figure |7| 
the horizontal dotted line denote the critical index n cr i t , corresponding to the point at which dX/dD = for a given A in 
(A, D)-plane. According to the prediction from the non-extensive thermostatistics with Tsallis entropy, after reaching the 
critical index n cr it, the system is expected to enter the gravothermally unstable regime. That is, the relaxation timescales 
between the inner and the outer parts are decoupled and the system finally undergoes the core-collapse. In a rigorous sense, the 
thermodynamic prediction itself cannot be applicable to the out-of-equilibrium state, however, it turns out that the predicted 
value n cr it or D C rit is consistent with the TV-body results. In fact, the successful fit is obtained until t ~ 34i r h,i, while the 
core-collapse takes place lately at t ~ 44i r h,i (see FiglSfr. Thus, the prediction based on the non-extensive thermostatistics 
may provide a crude estimate of the boundary between the stability and the instability in the general non-isothermal states. 
This point will be further discussed in other runs. 



4-2.2 Other runs 

Figures ISl ll II show the snapshots of the density profile, the distribution function and the velocity dispersion profile obtained 
from the runs n5A-n5D, which were all started with the same polytrope index n — 5, but with the different initial density 
contrasts; D = 100(run n5A, FigHJ: D = 10 3 (run n5B, FigEJ, D = 10 4 (run n5C, FigEHJ and D = 10 6 (run n5D, FiglTTl. 
Also in figures 1131 and 1141 the results taken from the runs nQA and n6B are depicted, whose initial density contrasts are 
D — 110 and 10 4 , respectively. 

Similarly to the run n3A, the results from the runs n5A-n5C , n6A and n6B exhibit the quasi-equilibrium regime in the 
early phase of the evolution. The evolutionary sequences in the quasi-equilibrium regime can be quantitatively characterized 
by a family of stellar polytropes with time-dependent polytrope index. Figures WI\ and [T5l summarize the fitting results for 
polytrope index as function of time in logarithmic scales. It seems apparently that for the cases starting with smaller value 
of D (i.e., runs n5A and n6A), time variation of polytrope index is systematically large and statistical fluctuation become 
noticeable, however, it turns out that this fact simply comes from the geometrical reason for the equilibrium sequence plotted 
in figure H That is, in (A, D)-plane, a number of trajectories of the stellar polytropes with different polytrope index n are 
assembled at the narrow region with smaller value, D or A, causing a large variation and/or fluctuation in the time evolution 
of polytrope index. 

Apart from this behavior, figures 1121 and 1151 show that the timescales of quasi-equilibrium state crucially depend on 
the initial conditions, which are basically characterized by the initial density contrast D or the dimensionless energy A. 
Qualitatively, the timescale of quasi-equilibrium state is understood from the local estim ates of the relaxation time, which are 
inversely proportional to the local density fe.g.. lSpitzerlll98^ : lBinnev fc Tremainelll987h : 

v 3 

*r = 0.065 — — . (29) 

G 2 mp In A 

For a system with small initial density contrast, the depth of the potential becomes shallower and equation 12911 implies that 
the relaxation proceeds slowly enough in both core and halo. The quasi-equilibrium state is thus expected to be long-lived. 
As anticipated in figures 151 and H31 fitting to the stellar polytropes is successful until t — 63f r h,i for run n5A and 66i r h.i for 
run n6A, much longer than the other cases. 

On the other hand, for a system with large initial density contrast, equation 1291 suggest that the decoupling of the 
timescales between core and halo would occur earlier and the quasi-equilibrium phase would be short-lived. Indeed, the 
distribution functions for runs n5C and nQB show that while the high-energy part of the function /(e) almost remains the 
same, the low-energy part of the distribution is rapidly developed and stretches toward e — > — oo. Also, the inner part of 
the density profile grows rapidly and shows an intermittent behavior with large fluctuation. As a result, fitting to the stellar 
polytropes failed earlier at t ~ 20t r h,i for run n5C and 12t r h,i for run n6B. After that, the system soon becomes unstable 
state, finally undergoing the core-collapse, consistent with the thermodynamic prediction based on the Tsallis entropy. Note 
that the fitting results in the runs n5C and n6B are slightly worse. From bottom panels of figures WI\ and [T51 the estimated 
X 2 values given by 128H become 5 ~ 14 in both runs, which is larger than the typical value 3 — 7 in cases with the smaller 
initial density contrast. 

The observation in both fitting results and the timescales for quasi-equilibrium states indicates that for some large values 
of the initial density contrast, the fitting to stellar polytropes would fail from the beginning and the quasi-equilibrium state 
ceases to exist. Indeed, such an example was obtained from the run n5D (see Fig lllH . The initial density contrast of this 
run is D = 10 6 and the dimensionless energy is A = 2.345. That is, the location of the initial state in the (A, D)-plane is 
outside the region depicted in figure H In this case, the central part of the system is highly concentrated and the core-halo 
structure is developed from the beginning. The resultant relaxation timescale is much shorter than that of the other cases 
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Figure 5. Lagrangian radii as function of time obtained from the run nSA. Each line represent the radius of the shell which contains a 
constant mass fraction. The numerical values of the mass fraction are indicated in the right part of figure. The horizontal dotted line is 
the potential softening scale, e = l/N ~ 4.9 X 10 — 4 . 
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Figure 6. Snapshots of density profile p(r) (Left), distribution function /(e) as function of specific energy of particles e = v 2 jl + $(r) 
(Middle) and one-dimensional velocity dispersion a 2 1D ()*) (Right) from the run n3A. In each panel, symbols represent the TV-body 
results, while the solid lines show the Emden solutions of initial distribution. The lines except for the solid lines are the results fitted to 
the one-parameter sequence of stellar polytropes. Note that for clarify, the results are offset vertically, successively by two-digits below, 
except for the final output. 
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Figure 7. Time evolution of polytrope index fitted to the ./V-body simulations in the case of the run n3A. The upper panel shows 
the fitting results for the runs with particles N = 2,048(solid), 4, 096(long- dashed) and 8, 192(short- dashed). The horizontal dotted line 
represents the critical index n cr ; t , corresponding to the marginal state of the entropy S q satisfying the condition dX/dD = in the case 
of run n3A. The lower panel represents the goodness of fit by plotting the function x 2 defined in equation 1281 . 
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and the system soon becomes unstable, leading to the earlier core-collapse. Compared to the stable stellar polytropes with 
the same initial energy A = 2.345 and with D < D C rit, none of the model parameters successfully reproduce the simulated 
density and/or velocity dispersion profiles (long-dashed, short-dashed and dot-dashed lines in Fie llf \ . Though the collapse 
time crucially depends on the softening parameter, the convergence test as examined in appendix A suggests that the collapse 
time converges to 18i r h,i, close to the standard result without the adiabatic wall, i.e., tcoU — 16t rh ,i (e.g., Tablel3.2 of lAarsethl 
120031) . 

Therefore, for a general initial condition with large D or A, quasi-equilibrium behavior generally ceases to exist. In other 
words, long-lived quasi-equilibrium states appear at the smaller value of A or D and the system is quantitatively characterized 
by the stellar polytropes. Their lifetime is expected to become much longer as approaching down to A = 0.335 and/or D = 709, 
i.e., the critical values for the marginal stability in isothermal distribution. 

4-2.3 Discussion 

So far, we have focused on the characterization of the quasi-equilibrium state using the one-parameter family of the stellar 
polytropes. While most of the transient state is well-approximated by the stellar polytropes with varying polytrope index, one 
may criticize that the use of the stellar polytropes is not the best characterization for the quasi-equilibrium sequence. Indeed, 
even restricting the stellar distribution to the stationary solutions of the Vlasov equation, one can, in principle, construct the 
infinite set of one-parameter family of stellar models. In this sense, the stellar polytropes should be regarded as a particular 
set of stellar models. Further, there is no rigorous proof for the uniqueness to characterize the quasi-equilibrium states. 

In real astronomical systems such as globular clusters, the stellar polytropes had been known to poorly fit to the observed 
structure of stellar distribution. Rather, the majority of the Galactic globular clusters is quantitatively characterized by the 
King model jKindll966t iBinnev fc Tremaindll987l : ISpitzerlll987l : iMevlan fe Heggielll997h . In contrast to the stellar polytropes 
as q-exponential distribution, King model is represented by the truncated exponential distribution: 

m oc { f - 1 ; < J , m 

where we define e' — e — </>(rs) with radius tb being the truncation radius and e is the specific energy of a particle. Provided 
the dimensionless energy A, the equilibrium sequence of the King model is then characterized by the one-parameter Wo = 
2/3[(/>(rs) — 0(0)], which represents the depth of the gravitational potential. 

In figure 1161 the TV-body data taken from the run n3A is used to compare with the King model. Note that the fitted 
values of the parameter Wo indicated by figure [TrU were obtained under the suitable restriction, rs > r e . Similar to the 
stellar polytropes, the density and the velocity dispersion profiles reasonably fit to the King model. The fitted value of the 
parameter Wo gradually increases as time goes on, indicating that the depth of the poten tial becomes deeper. Compared with 
the observed Galactic globular clusters with typical range Wo = 4-10 llTrager et al.lll995l) . fitting results for W are somewhat 
large. 

On the other hand, turn to focus on the distribution function, deviation from the King model becomes manifest at the 
high-energy tails e > 0. While the distribution function for King model falls off at the truncation energy e = ^{tb) which takes 
the negative value, number of high-energy particles with e > gradually increases in our TV-body calculation, which contributes 
to a high-energy tail of the distribution function. Although the low-energy part of the distribution function resembles the 
exponential form of the King model, the discrepancy at the high-energy part implies that the boundary condition in our 
idealistic situation is very different from that of the Galactic globular clusters. In fact, the influence of external tidal field 
is significant for the Galac tic globular cluster and the resultant distribution function sharply falls off at the tidal boundary 
jBinnev fc Merrifieldlll998ft. The stellar particles which are usually bounded to the system tend to esc ape from the globular 
cluster system (e.g.. lFukushige fc Heggiell995l : lBaumgardt fc Makinol2003l:lTanikawa fc Fukushigel2005l) . To mimic this effect, 
truncation radius vb is artificially introduced by hand in the King model. On the other hand, in presence of the adiabatic 
wall, the high-energy particles with e > 0, which are usually unbounded, cannot freely escape outward from the system and 
thereby no specific energy cutoff appears. 

As a result, the stellar polytropic distribution which has no energy cutoff successfully reproduces the quasi-equilibrium 
states in our TV-body setup. This means that, in presence of the adiabatic wall, the simple power-law distribution provides a 
better characterization than the truncated exponential distribution. Since the adiabatic wall is an artificial but the simplest 
boundary condition, the stellar polytropic distribution might be regarded as a fundamental stellar model theoretically, though 
not practically useful in characterizing the observed structure of Galactic globular clusters. 

4.3 A family of stellar models with cusped density profile 

The analyses in previous subsection have revealed that the quasi-equilibrium evolution characterized by the stellar polytropic 
distribution can appear at the energy A close to the critical value of the isothermal distribution. We then attempt to clarify the 
generality and/or the physical conditions for the quasi-equilibrium state in more general initial conditions. In this subsection, 
as a special cla ss of initial conditions that contain the non power-law features, we treat a family of stellar models with cusped 
density profile llTremaine et al.lll994h . The models contain two parameters, one is the scale-radius a, and the other is related 
to the slope of the inner density profile r\. The density profiles of these models are expressed as: 
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Figure 9. Same as Fig|S] but in case of the run nbB. 
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Figure 10. Same as Fig|g] but in case of the run n5C. 




Figure 11. Same as Fig[5] but in the case of the run n5D. 
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P( r ) K / , SI / Mi ■ (31) 

Note that the above expression includes the models considered bv lHernauistl il990l) for r\ = 2 and bv I Jaffd 1^983) for rj = 1 
as special cases. 



4-3.1 N -body setup and overview 

For a spherically symmetric configuration with isotropic velocity distribution, usi ng the Eddington formula , the one-particle 
distrib ution function can be uniquely reconstructed from the density profiles (Bin nev fc Tremaindll98^ . lTremaine et alJ 
l)l994F) gave useful analytic formulas for various physical quantities such as the distribution function, the p otential and the 
velocity dispersion profiles for a system extending over the infinite radius, r — > oo. Based on their formulas. iQuinlanl |l996) 
numerically estimated the dependence of the slope r\ on the timescale of core-collapse using the Fokker-Planck code. For 
present purpose, however, a direct application of their formulas is inadequate in presence of the adiabatic wall. 

In appendix B, taking account of the truncation radius r e , we re-derive the analytic formulas for distribution functions 
as well as the other physical quantities. Based on this, figure [T7I plots the theoretical curves for the distribution function(Ze/t) 
and the velocity distribution profile (right) with a specific choice of the parameters, r\ = 1, 1.5, 2 and 3, fixing the scale-radius 
to a/r e = 0.5. For models with r\ > 1, the distribution function exhibits a divergent behavior at a finite energy e = e m i n < 0, 
/(emin) = +oo. Further, the velocity dispersion profile shows non-monotonic behavior; it first increases and eventually turns 
to decrease as approaching the center. These peculiar features imply that the inward heat flow occurs along a course of the 
relaxation, causing the central part of the system expand, which is the same phenomenon as observed in the unstable case of 
the isothermal distribution (see Sec 14. in . 

To perform a simulation, we increased the number of particles to N = 8K in order to resolve the central part of the 
cuspy density profiles. The softening parameter of gravitational potential is set to e = A/N. Note that the convergence test in 
appendix A suggests that much smaller value of the softening parameter should be used for a system with highly concentrated 
core. However, decreasing the softening parameter requires a much longer calculation time and the probability of binary- 
formation via three-body interaction increases around the core. Since we use a simple individual time-step algorithm without 
any regularization schemes, the TV-body integration becomes heavily time-consuming once a tight binary is formed. Hence, we 
do not discuss here the timescale of the quasi-equilibrium state and focus only on the condition for quasi-equilibrium states. 
The quantitative estimates of the timescale will be presented in future task. 

In table |3] the model parameters of the initial conditions examined here are summarized, together with the evolutionary 
status. Also in figure tTHl the result for each run is represented by the symbols in energy-scale radius relation. Roughly speaking, 
figure^] says that the quasi-equilibrium sequence characterized by the stellar polytrope appears(open stars, case(B)) when 
the total energy of the system is not far from the critical energy of isothermal distribution, A cr i t = 0.335. However, this is 
not a sufficient criterion for the presence of quasi-equilibrium state. For models with r\ = 1, whose density profile resembles a 
singular isothermal sphere at the center, the system cannot be fitted by the stellar polytrope (filled stars, case(A)). On the 
other hand, for the initial conditions with A < A cr i t (stars in shaded-region) , the system finally approaches the isothermal 
state. In this case, the transient state of the system could be also fitted by the stellar polytropes, although the fitted value of 
the polytrope index is so large that one cannot be easily discriminate it from the isothermal distribution. 
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Figure 13. Same as FiglHI but m the case of the run n%A. 






Figure 14. Same as Fig[5] but in the case of the run n6B. 
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Figure 15. Evolution of polytrope index in the cases of runs n6A and n6B 



4-3.2 quasi- attractive behaviors and condition for quasi- equilibrium state 

Let us focus on the characteristic behaviors of the long-term evolution by picking up some typical examples. Figures 1191 1201 
and !21l show the results obtained from the runs r)2C, r/1.5B and rjlC. Figure EH shows the typical example in which the system 
finally approaches the stable isothermal state. Due to the small value of A, the effect of self-gravity is small and the system 
evolves very slowly. While we tried to fit the transient states of the system to the stellar polytropes, the resultant value of 
the polytrope index n is quite large, n ~ 20-60, indicating the system being in nearly isothermal state. For comparison, we 
also plot the theoretical curves for isothermal distribution. It seems difficult to discriminate which models are fitted to the 
simulation results better. 

By contrast, in figure I2T1 the quasi-equilibrium state approximated by the stellar polytropes appears. In the run r/1.5B, 
because of the dimensionless energy A > 0.335, the system finally undergoes core-collapse. Looking at an early phase, however, 
the core expansion first takes place and the flatter core is formed. Then the core density turns to increase gradually and the 
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Figure 16. Same as Fig|5| but results are compared with the King models. The fitting results shown in the lines except for solid lines 
are obtained by varying the dimcnsionless parameter, Wo = 2f3[(p(rg) — <j>(0)] under keeping the mass and the total energy fixed. 

transient state can be approximately described by the stellar polytropes for a long time (t ^ 30ttr,i). These behaviors, which 
may be regarded as quasi-attractive behavior, can be explained from the inner structure of the initial velocity dispersion 
profile. That is, due to the small value of the local relaxation time t r at the core, the inward heat flow first occurs toward the 
equipartition of the kinetic energy, leading to the uniform velocity dispersion at the inner part. Then, the outward heat flow 
next occurs and the inhomogeneity in the velocity dispersion is gradually erased. Although this slow relaxation does not stop 
and finally leads to the catastrophic heat flow, the system is remarkably long-lived. 

On the other hand, the run r/lC has slightly smaller value of A than the run r/1.5B and one naively expects that the 
system is stable. However, the results shown in figure I2T1 is completely opposite. In contrast to the run T/1.5B, the initial 
condition of the run r/lC has uniform velocity dispersion at the inner part. This implies that the inward heat flow is only 
supplied by the randomness of the initial perturbation, as seen in the unstable isothermal case (see Sec l4.lt . Thus, the amount 
of the heat flow is insufficient and the resultant core radius is rather small, whose density profile cannot be approximated 
by the stable stellar polytropic distribution. For comparison, in figure I2T1 we plot the stellar polytrope with index n ~ 17.6, 
which is the marginal stable state that has the same total energy A as in the run r/lC. While the distribution function and 
the velocity dispersion profile resemble the marginal stable state of the stellar polytrope, the discrepancy is manifest in the 
density profile. As a result, the system is short-lived and could not reach the quasi-equilibrium state. 

Note that the lifetime in the model nlC sensitively depends on the randomness of the initial perturbation. Figure 12 1 1 
shows the run-by-r un variation of the time evolution of the core radii, where th e core radius was estimated according to the 
procedure given bv lCasertano fc Hud (ll985Tl (see also iFukushige fc Heggielll995P . Compared to the runs rflC and nl.5B, the 
run-by-run variation in nlC is significant and it seems to originate from the first stage of the core expansion. This behavior also 
holds for the other runs nlA and nlB. Therefore, the lifetime of the system starting from the initial conditions with cusped 
density profile p oc r~ 2 would be generally stochastic. Although the present calculation with non-zero softening parameter 
is inadequate to estimate the precise core-collapse time, the uncertainty of the collapse time would remain true even if the 
appropriate regularization scheme is implemented in our code. In other words, the condition for quasi-equilibrium state as 
well as the lifetime of the system are much sensitive to the velocity structure of the initial condition. Though the present 
surveys do not give a conclusive statement for generality of quasi-equilibrium state, the out-of-equilibrium state with the 
cusped profile p oc r~ a , (a < 2) possibly exhibits the quasi-equilibrium behaviors if the dimensionless energy A is close to 
0.335. 



4-3.3 Entropy growth and quasi-equilibrium state 

In order to elucidate the sequence of the quasi-equilibrium evolution from a thermodynamic point-of-view, we quantify the 
entropy growth for the A-body data, r]2C, r/1.5B and r\\C . Though the quasi-equilibrium behavior seen in the simulations 
may imply that the entropy growth can be better characterized by the non-extensive Tsallis entropy, the entropy measure in 
quantifying S q is not, strictly speaking, unique when comparing with different values of the polytrope indices n, or equivalently, 
the q-parameters. Therefore, as long as the cases with time-varying polytrope indices are concerned, it would be natural to 
quantify the entropy growth with the Boltzmann-Gibbs entropy. 

In left panel of figure |2"31 the results are plotted as the trajectories in (D, 5 , B G/A)-plane, together with the equilibrium 
sequence for the stellar polytropes denoted by continuous lines. The derivation of theoretical curves for isothermal and the 
stellar polytropes is described in appendix C. Also, in right panel, the entropy growth is quantified and is plotted as function 
of time. Note that the time interval between the symbols marked along each trajectory roughly corresponds to a half-mass 
relaxation time for the initial distribution. 

In figure 12^1 while the time evolution of density contrast D shows non-monotonic behaviors, the specific entropy Sbg/A 
monotonically increases in time, consistent with the law of thermodynamics indicated by the Boltzmann H-theorem. The 
evolutionary sequence in the (D, S l BG/A)-plane depends on the initial condition. In the case of the run r]2C, the left panel 
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Figure 17. Distribution function(Le/4) and velocity dispersion profile (Right) for a family of stellar model by Tremaine et al. (1994) in 
presence of an adiabatic wall. For a specific value of the scale radius a/r e = 0.5, the analytic results with various slopes r\ are plotted 
based on the formulae in Appendix B. 



Table 3. Model parameters and the evolutionary states in cases starting from the stellar models bv lTremaine et alj ll99 



run # 




parameters 


half-mass radius (rh/r e ) 


# of particles 


realizations 


transient state 


final state 


T)1A 


V 


= 1, a/r e = 0.5 


0.25 


8k 


2 


none 


collapse 


r/lB 


V 


= 1, a/r e = 0.8 


0.3077 


8k 


6 


none 


collapse 


r/lC 


V 


= 1, a/r e = 2.0 


0.4 


8k 


4 


none 


collapse 


771.5A 


V = 


= 1.5, a/r e = 0.2 


0.221 


8k 


2 


none 


collapse 


771.5B 


V = 


= 1.5, a/r e = 0.5 


0.362 


8k 


3 


stellar polytrope 


collapse 


771.5C 


V = 


= 1.5, a/r e = 1.0 


0.4598 


8k 


3 


stellar polytrope 


collapse 


r)2A 


1 


= 2, a/r e = 0.2 


0.2869 


8k 


2 


none 


collapse 


772B 1 


n 


= 2, a/r e = 0.5 


0.4459 


8k 


2 


stellar polytrope 


collapse 


r/2C 2 


V 


= 2, a/r e = 1.0 


0.5469 


8k 


3 


none 


isothermal 


r/3A 


1 


= 3, a/r e = 0.2 


0.3907 


8k 


2 


stellar polytrope 


collapse 


r/3B 


1 


= 3, a/r e = 0.5 


0.5619 


8k 


2 


none 


isothermal 



1 initial distribution corresponding to run CI in 

2 initial distribution corresponding to run C2 in 



Taruva fc Sakagani] <2003d 
Taruva & Sakagami 1 2003c 



shows that the trajectory is located around Sbg — 3.7 — 3.9 and the evolution slows down after contacting with the trajectory 
of the isothermal distribution. As for the trajectory of the run r/1.5B, it first goes across the boundary between the stable and 
the unstable stellar polytropes indicated by the dotted line. Then, it temporarily settles down into a stable stellar polytropic 
state with index n ~ 12 — 14 for a long time. At that time, the transient states were successfully fitted by the stable stellar 
polytropic distribution (see Fig |2UH . It is interesting to note that the growth of the entropy shown in the right panel of 
Figure 12^1 is slightly restrained during the quasi-equilibrium regime, which might manifest the minimum entropy production 
principle in non-equilibrium thermodynamics. On the other hand, for the trajectory of the run r]lC, while it first approaches 
the stability boundary, the decreases of the density contrast eventually terminate at relatively higher value D ~ 10 4 and 
accordingly the transients could not be fitted by the stable stellar polytropes. 

Note also that the transient state of the run r/lC cannot be even fitted by the unstab le stellar polytropes, since the 
density profile of the unstable polytropes shows a log-periodic behavior at the outer part (e.g.. IChand rasekhar 1939) while no 
such behavior appears in the A^-body simulation. In this sense, only the quasi-equilibrium state characterized by the stable 
stellar polytropes might have some special meanings. Similar to the isothermal distribution, the quasi-equilibrium state has 
the quasi-attractive feature that the system starting from some classes of initial conditions tends to approach the polytropic 
state, which may provide an important suggestion for the reality of the non-extensive statistics. 



5 DISCUSSION & CONCLUSION 

In this paper, we have discussed a possible application of non-extensive thermostatistics to the stellar systems and attempted 
to explore its reality. To do this, we have numerically investigated the quasi-equilibrium properties of the A-body systems 
before the core-collapse stage. Particularly focusing on the long-term stellar dynamical evolution from the thermostatistical 
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Figure 18. Energy-scale radius relationship in stellar models by Tremaine et al. (1994) in presence of an adiabatic wall. The thin lines 
represent the equilibrium sequences of stellar model with different values of index r). Along each line, we plot the star symbols, which 
indicate the result of Af-body simulations fitted to the stellar polytropes. The fitting results shown in the figure are categorized as case 
(A) and case (B), which are roughly separated by thick dashed line. The case (A) denoted by filled symbols means that fitting to the 
(stable) stellar polytropes was failed, while the case (B) denoted by open symbols indicates that the transient states approximated by 
the stellar polytropic distribution appeared as quasi-equilibrium state. 

point-of-view, we try to characterize the out-of-equilibrium state starting with various initial conditions in the setup of the 
so-called Antonov problem. We found that the quasi-equilibrium states, in which the system evolves gradually on timescales 
of two-body relaxation, appears and the system may follow the quasi-equilibrium sequence for a long time if the dimensionless 
energy A = —Er e /GM 2 of the system is not so far from the critical energy of the isothermal sphere, A cr it = 0.335. The 
schematic illustration of our basic results is shown in figure I2H The transient states during the quasi-equilibrium evolution are 
approximately described by the one-parameter family of stellar polytropes as extremum states of non-extensive entropy S q with 
the time- varying polytrope index. The fitted value of the index n gradually increases with time and the system keeps following 
a sequence of stellar polytropes until reaching the critical index, n cr it- In general, the condition for the quasi-equilibrium state 
would depend on the details of the velocity structure in the initial conditions, however, within a class of initial conditions 
examined in this paper (i.e., stellar polytropes and stellar models with cusped density profiles), the out-of-equilibrium states 
with inner density profiles p oc r~ a (a < 2) (or rj ^ 1) exhibit the quasi-equilibrium behavior that is attracted to a sequence 
of stellar polytropes. 

One may naively think that the results obtained here severely depend on the presence of an adiabatic wall, since the real 
stellar systems in absence of the adiabatic wall are known to be poorly fitted to the stellar polytropes. Recalling the discussion 
in section 14.2.31 however, the outer part of the system is expected to be mainly affected by the boundary condition, since the 
relaxation timescale at the outer part is rather longer than that at the core. In other words, as long as the relaxation time 
at the central part is shorter than that at the outer part, the modification of the boundary condition only alters the outer 
part of the system, not all of the system. In fact, we have seen in section 14.2.31 that the phenomenological King model which 
accounts for the globular clusters affected by the Galactic tidal field resembles the stellar polytropes at the inner part. In 
this sense, the stellar polytropic system as quasi -equilibrium state would be a fundamental stellar model and may sometimes 
make sense even if removing the adiabatic wall jTaruva fc Sakagamilteoolkl) . 

The present results may give an interesting suggestion for the justification and/or the realization of the non-extensive 
thermostatistics based on the Tsallis entropy. Strictly speaking, however, the A^-body simulations just indicate a reality of the 
stellar polytropes as q-analogue of the Boltzmann-Gibbs dist ribution. Further, the q uasi-equilibrium state is time-dependent, 
which cannot be rigorously treated by the thermostatistics iChavanis fc Siral2004l) . Even at this moment, it might be also 
possible to give the interpretation that stel lar polytrope is merely a dynamical equilibrium state, but it can be stable even 
with respect to a non-linear perturbation iChavani s 2003). In this respect, the present A^-body results are not conclusive 
and the interpretation of their results should be carefully discussed. Nevertheless, one may hope that the exploration of a 
connection with non-extensive entropy opens a new window to understand the non-equilibrium thermodynamics of the long- 
range systems. In this respect, the analytical treatment based on the kinetic theory is an important next step to interpret 
the A^-body results thermodynamically. A crucial task is to estimate the timescales for quasi-equilibrium evolution, as well as 
to determine the generic criteria for quasi-equilibria. To investigate this, the Fokker-Planck model for stellar dynamics would 
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Figure 19. Snapshots of density profile, distribution function and velocity dispersion for run r]2C. 




Figure 20. Snapshots of density profile, distribution function and velocity dispersion for run r)l.5B. 




Figure 21. Snapshots of density profile, distribution function and velocity dispersion for run r\\C. 



be helpful iTaruva fc Sakagamil (l2004ll . A quantitative comparison between a numerical solution of the Fokker-Planck model 
with the iV-body results based on a more sophisticated iV-body code such as the one developed by Aarseth will be presented 
elsewhere. 
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Figure 22. Run-by-Run variation of evolution of core radius among three different realization for the runs nl.5B and r]2C, and among 
four different realization for the run r/lC. The dotted line indicates the length of the potential softening, e = 4/N ~ 4.9 X 10 — 4 . 




Figure 23. left: Trajectories in (SbGi D )-plane obtained from the runs rjlC (open- squares), r]1.5B(filled triangles) and rj2C (open- circles). 
The lines with symbols indicate the simulation results, while the solid, the short-dashed, the long-dashed and the dot-dashed lines are the 
equilibrium sequences of stellar polytropes with different polytrope index n, which are obtained from the analytic results in Appendix C. 
On the other hand, the dotted curve represents the stability boundary of the stellar polytropes, inferred from the energy-density contrast 
relation (dotted curve in Fig|5J- Right: Time evolution of Boltzmann-Gibbs entropy from the runs nlC, ijl.5B and r]2C. 



APPENDIX A: CONVERGENCE TEST OF N-BODY SIMULATION 

In a series of our iV-body simulations, the Plummer softening with parameter e = 1/N or 4/N is used to avoid the formation 
of tight binaries. For our investigation of the quasi-equilibrium state before the core-collapse stage, the fourth order Hermite 
integration code with individual time-step provides a robust numerical method without a regularization scheme. However, it 
would be crucial to pursue the gravothermally unstable regime. While our primary concern is the non-equilibrium evolution 
before entering the core-collapse phase, it is important to note the effect of potential softening on the estimation of collapse 
time. 

To quantify this, convergence property of the collapse time is investigated by varying the softening parameter e. For this 
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Figure 24. Schematic illustration of the simulation results, which summarizes the evolutionary status shown in figures l2l and Il8l 



purpose, TV-body simulations starting with the Plummer model, i.e., the stellar polytrope with index n — 5, are used to study 
the effect of potential softening. The corresponding simulations are the runs n5B and n5D listed in table |2] 

Figure Q shows the results of the convergence test based on the run n5B. Left panel plots the evolution of core radii for 
various choices of softening parameter, while the right panel estimate the collapse time as function of softening parameter 
e for three different realizations of the initial condition. Here, the collapse time is defined as the first passage time when 
the core radius becomes shorter than 3e. Clearly, large value of the softening parameter overestimates the estimation of the 
core-collapse time and the transition between the quasi-equilibrium stage to the core-collapse phase becomes uncertain. For 
an appropriately small value e £ 1/(AN) ~ 1.2 x 10 -4 , the core radius sharply falls off around the core-collapse time and the 
collapse time seems to converge to t co n ~ 30-40t r h,i, although there exists a large scatter among three different realizations of 
the initial condition. 

In general, the convergence properties of the collapse time depend on the initial conditions. Figure shows the results 
from the run n5D, i.e., Plummer model with density contrast D — 10 6 . In this case, the scatter becomes even smaller and the 
collapse time converges to t co ii ~ 18£ r h,i at e ^ 6 x 10~ 5 , slower than the case of run n5B. These experiments indicate that 
a large softening length affects not only the core-collapse time but also the timescales of quasi-equilibrium stage. Further, 
the requirement for the softening parameter becomes severe as increasing the initial density contrast. Hence, for quantitative 
estimation of collapse time, the softening length should be set to zero with implementing a regularization scheme. The use of 
a small softening length is also favorable to investigate the quasi-equilibrium evolution. FiguresQand|2]suggest that for initial 
conditions with moderate range of the initial density contrast, 10 2 ^ D ^ 10 4 , the softening length with e ^ l/N provides a 
better estimation of the quasi-equilibrium timescales, although it still overestimates the collapse time. Hence, in this paper, 
we mainly use the softening length e = l/N. 



APPENDIX B: A FAMILY OF STELLAR MODEL WITH CUSPS 

In this appendix, we present analytic formulae for the stellar model considered bv lTremaine et all (1994) taking account of 
the adiabatic wall. 

Just for convenience, let us first introduce the following variables: 

M n GM fr e + ay 



Alia 3 ' 

In terms of these, the density profile for stellar models bv lTremaine et all ill 9941) is expressed as 

^ = (r/a) 3 -i(l + r/a) 1 +*'' ^ 
The corresponding mass profile and the gravitational potential respectively becomes 

m(r) = / dr' An r' 2 p(r') — — ( — - — \ , (3) 
Jo V \r + aj 




Figure 1. Convergence properties of collapse time for various choices of softening parameter e in the case of run n5B. Here, the collapse 
time is defined as the first passage time when the core radius becomes shorter than 3e. Left: evolution of core radius for a realization with 
seed number 327. Right: Variation of collapse time as function of softening parameter among three realizations with different random 
seed numbers. 




and 

$(r) = -G { + / dr Anr'p(r') 



-B 



-B 



-lVr + J + (r e /a + lW ' V ^ 



re/a)" 

rj — l\r e + aj rj— l\r + aj ' (r e /a + 1) 1 ' 

(4) 



{log(^)-log(^) + -i— 1 ; , = 1 

\r e + aj \r + aj r e /a+l) 



For stationary state of the Vlasov equation, the one-dimensional velocity dispersion profile can be calculated from the 
density profiles and the gravitational potential through the Jeans equation. In the case of the isotropic velocity distribution, 
this gives 



< 1D (r) = (5) 
with the function P{r) being pressure, determined from the hydrostatic equilibrium(Jeans equation), dP/dr = — p d<&/dr. 
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Note that the explicit expression for function P(r) cannot be expressed in an unified manner. For some values of 77, we have 



P„=i.5(r) 



AB 
AB 
AB 
AB 



6 log 
4 log 
4 log 



V + l 

y 
y 

y + i 

y + i 
y 



+ 



6(l + y)-9 + 
2 



y+1 2(y + l)2 
2 1 



y + 1 3(y + l) 2 3(2/ + 1)3 
11 1 



y+1 2(y + l) 2 3(y + l) 3 4(y + l)' 



(6) 
(7) 
(8) 
(9) 



_5(y + l) 5 6(y + 1) 

where we defined y = r/a. Provided the pressure, one can also calculate the total energy of the system confined in a wall as 



E = K + 



U = J d 6 d 3 xd 3 v ji?j 2 + f( r , v ) 



dr 4nr P(r) + — j dr Airr p(r)$(r). 



In terms of the dimensionless variable A = —r c E/(GM 2 ), we obtain 



A I)= i — 



^77=2 



- (1 + 5y e + 18y 2 + 12y 3 ) - 3y 2 (y e + l) 2 log 
~ (-3 - 43y e - 60y 2 - 24yf ) - 3(y e + l) 3 log 



Ve + 1 



= — (29 + 53y e + 42y 2 + 12y e 3 ) - (y e + l) 4 log 



Ve + l 
Ve+l 



Ve 



A 



17=3 = 



20y; 



2 (l + 7y e +J/ e 2 ) (tfe-1) 



(10) 

(11) 
(12) 

(13) 
(14) 



where the variable y e denotes r e /a. 

According to the standard text book for stellar dynamics, the one-particle distribution function for the isotropic spherical 
stellar model is expressed as a function of the specific energy e = v 2 /2 + $(r) and can be reconstructed from the density 
profiles through the Eddington formula. Introducing the variables e — $0 — e and ip = <f>o — ®(r), we have 



m = 



d p dtp 



V8tt 2 J dip 2 y/F^ip 
with the regularity condition: 



= 0. 



(15) 



(16) 



Note that the numerical constant $0 is determined from the above condition. With the use of the analytical expressions 
and @, after some manipulation, the Eddington formula 1151 can be rewritten with 



1 A 



9(e) 



dip 12{s(V>)} 2 ~ 4(t) -4)s(V>) + 2(3-r;) 



where the functions g(e) and s(V0 ar e respectively given by 



(r e /o) 



17-1 



q{e) = - 1 + ^-^ (r. - J i" J ' 



SW 1-{1- (77- l)^} 1 /^-!) 

in cases with 77 7^ 1 and 



(17) 

(18) 
(19) 

(20) 
(21) 



s(ip) = — 

for 77 = 1. In principle, the distribution function /(e) is obtained from the numerical integration of 1151 . In some specific 



24 A. Taruya and M. Sakagami 



values of 77, however, one can luckily obtain the analytic expressions of /(e): 



/r;=i.s(e) 
A=s(e) 



.4 



5 3 / 2 V87T 2 (2-<?) 9 /2 



-(3 + 32g-8g 2 )sin _1 



Vg(2 - g) 

28 

A 1 



{63 + 693g - 5670g 2 + 7410g 3 - 4488g 4 + 1448g 5 - 240g 6 + 16g 7 } 
1 



5 3 / 2 4^2 (1 - g) 5 / 2 
4 1 



3 HhT^Vg) - V^O^g) (16g 3 - 24g 2 + 2g + 3) 



2A/2g (3-4g)+3(l-2g) log 



53/2 n 2(! _ 2q ) 

Here, the function F(x) and erf(a;) are Dawson's integral and the error function: 



erf (a;) = V2-K / die 



F{x) = e 



t 

dt e 



(22) 



(23) 
(24) 
(25) 

(26) 



Compared the final expression s 12211 - 1251 with those obtained bv lTremaine et, all l)l994f) . we deduce that the only alternation 
in the expressions of lTremaine et all (^994) in presence of the adiabatic boundary is to replace all the variables e at the right- 
hand side of equations with the function g(e) defined above. Therefore, the above results consistently recover the formulas 
derived bv lTremaine et all (|l994j) in the limit r e — *• 00. 



APPENDIX C: BOLTZMANN-GIBBS ENTROPY FOR STELLAR POLYTROPES 

Here, we derive the explicit expression for Boltzmann-Gibbs entropy in the case of the stellar poly tropic distribution, which 
is shown in theoretical curves of figure l2*3l 

First, we substitute the stellar polytropic distribution [1 1 lit into the definition of Boltzmann-Gibbs entropy: 

s^/n = - d 3 xd\ (L—)m' J 



N J \ N 

= -(sa + sb + sc)- (27) 
The quantities sa, sb and sc involving the phase-space integral are separately evaluated as 

sa = —= In I —= : \ I d*xd s v Pth !, „,„ (l - v2/2 " 



4^2tt 5(3/2, n- 1/2) |_ V2tt 5(3/2, n - 1/2) J J ' {(n + l)5/p} 3 / 2 [ (n + l)P/p 

1 

4V2nB(3/2,n - 1/2) 

sb = — -= — — / d 3 xd 3 v — ^— —J- I 1 - 1— — \ In <M 



4^5(3/2,71-1/2)7 {(n + l)P/p} s / 2 \ (n + l)P/p J \ (n + l)P/p 

3 



= (n-|){^(n-l/2)-^(n + l)} I 
= (n-|){V(n-l/2)-V(n + l)}, (29) 

sc = I /"^ L v 2 /z \ n - 3/ \ P/E 

4>/27r 5(3/2, n- 1/2) J ' {(n+l)P/p}W \ (n + l)P/pj {{n + l)P/p}^ 

= [fxJLin — Pi** (30) 

J M {(n + l)P/p} s / 2 v ' 

Here, the functions ip{z) and B(a,b) are the digamma and the beta functions, respectively. Collecting the above terms, the 
Boltzmann-Gibbs entropy becomes 

Sr >> = B (§, n - ») , - („ - ') {* (. - i) - * (. + 1)} - / *. <L h ( T? _^L 



l)P/p}WS 

For the remaining integral in the above expression, a further reduction is possible by repeating the integration by part. 
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Using the equations of hydrostatic equilibrium, a straightforward calculation yields 

M ln {( n + l)P/p}3/2 - [nM \ 

n- 3/2 



(w + l)P e 

l+l/n 



In p e 



+ 



n + 1 P e r e 



M 



6 + 



(n + l) 2 



G 2 m 3 



Mr 3 



(31) 



Note that the subscript e denotes the quantities evaluated at the wall, r — r e . In terms of the homology invariants (u,v) 
defined in equations 12111 and (1221 . the entropy Sbg is finally expressed as: 



o(poiy) 

^BG 



/N = -6 



n - 3/2 



+ ln 



n 

n - 3/2 
n + 1 

|27r(n + l) 



(H) M^H^} 
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u e — In 



3/2 
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+ - ln(27rGMr e ) 



/2 B{-, 
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(?) 



Mr 3 



In the limit n — * oo, we correctly recover the Boltzmann-Gibbs entropy for the isothermal state: 



9 



+ v e + 2u e 



In 



3/2 
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+ -ln(27rGMr e ). 



(32) 



(33) 
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